A continuación se expone el código y comentarios resultados de la realización de la práctica 1.
library(tidyverse)
library(tidyr)
library(dplyr)
library(knitr)
library(ggplot2)
library(stringr)
library(reshape)
library(minerva)
library(heatmaply)
library(kableExtra)
library(factoextra)
library(ggbiplot)
library(randomForest)
library(ROCR)
library(ggcorrplot)
library(ModelMetrics)
Aunque se nos entrega el fichero de train entero y podriamos obtener la variable respuesta, se trataran cómo si de un caso real se tratara, y todo el EDA y el training se usara solo el p1_train. El p1_test se usara únicamente para validar los resultados al final de todo de la practica.
train <- read.csv("./data/train.csv", head = TRUE, stringsAsFactors=TRUE)
p1_train <- read.csv("./data/p1_train.csv", head = TRUE, sep=";", stringsAsFactors=TRUE)
p1_test <- read.csv("./data/p1_test.csv", sep=";", head = TRUE, stringsAsFactors=TRUE)
En este apartado se va hacer el análisi exploratorio de los datos.
head(p1_train)
Eliminamos la variable id del dataset ya que es el indice, y no aporta información.
p1_train <- p1_train %>% select(-id)
Summary de los datos
summary(p1_train)
year hour season holiday workingday weather temp atemp
Min. :2011 Min. : 0.00 Min. :1.000 Min. :0.000 Min. :0.0000 Min. :1.00 Min. : 0.82 Min. : 0.76
1st Qu.:2011 1st Qu.: 6.00 1st Qu.:2.000 1st Qu.:0.000 1st Qu.:0.0000 1st Qu.:1.00 1st Qu.:13.94 1st Qu.:16.66
Median :2011 Median :12.00 Median :3.000 Median :0.000 Median :1.0000 Median :1.00 Median :20.50 Median :24.24
Mean :2011 Mean :11.57 Mean :2.506 Mean :0.029 Mean :0.6773 Mean :1.41 Mean :20.27 Mean :23.70
3rd Qu.:2012 3rd Qu.:18.00 3rd Qu.:4.000 3rd Qu.:0.000 3rd Qu.:1.0000 3rd Qu.:2.00 3rd Qu.:26.24 3rd Qu.:31.06
Max. :2012 Max. :23.00 Max. :4.000 Max. :1.000 Max. :1.0000 Max. :3.00 Max. :41.00 Max. :45.45
humidity windspeed count
Min. : 0.00 Min. : 0.000 Min. : 1.0
1st Qu.: 46.00 1st Qu.: 7.002 1st Qu.: 41.0
Median : 62.00 Median :12.998 Median :145.0
Mean : 61.77 Mean :12.802 Mean :191.4
3rd Qu.: 77.00 3rd Qu.:16.998 3rd Qu.:283.0
Max. :100.00 Max. :56.997 Max. :977.0
Se puede apreciar que no encontramos ningún nan en las columnas del dataset, haciendo que no tengamos que tratar con ellos.
p1_train %>% group_by() %>% summarise_all(funs(sum(is.na(.))))
Se eliminan las filas duplicadas del dataset.
p1_train <- p1_train %>%
distinct()
p1_train
Para las variables que son categoricas, aunque su representación sea numerica debemos tratarlas cómo categoricas a nivel dato para obtener el resultado deseado.
Las variables categoricas del dataset segun el enunciado son: workingday, holiday, weather y season.
p1_train <- p1_train %>% mutate(workingday = as.factor(workingday),
holiday = as.factor(holiday),
weather = as.factor(weather),
season = as.factor(season))
Eliminamos la variable target del dataset para hacer un estudio de las columnas predictivas. Posteriormente se hará el estudio de la variable target independientemente.
predictors <- p1_train %>% select(-count)
(Sobre variables numericas)
meltPred <- predictors %>% select(where(is.numeric)) %>% melt()
Using as id variables
meltPred %>%
ggplot(aes(factor(variable), value)) +
geom_violin(width=1, color = "gray", alpha = 0.2, fill = 'green' ) +
geom_boxplot(width=0.3, color="black", fill='blue', alpha=0.3, outlier.colour="red", outlier.shape=8,
outlier.size=1, notch=TRUE) + facet_wrap(~variable, scale="free")
Notch went outside hinges
ℹ Do you want `notch = FALSE`?
Estos gráficos representan un boxplot y un violinplot sobre todas las variables numericas del dataset. Podemos apreciar:
La columna year contiene únicamente dos valores 2011 y 2012.
La columna hour se ditribuye entre el 0 y el 24, con una media situada en los 12. No encontramos outliers.
La columna temp y la columna atemp siguen distribuciones muy similares, la media esta alrededor del 25. No encontramos outliers.
La columna humidity se distribuye entre el 0 y el 100 con una media alrededor de los 60. No encontramos otuliers.
La columna windspeed se distribuye entre el 0 y el 60 con una media alrededor de los 15. Encontramos numersos outliers cuando el valores es mayor de 30.
corrp <- p1_train %>% select(where(is.numeric)) %>% cor(method = 'pearson')
corrp %>% ggcorrplot(hc.order = TRUE, type = 'lower', lab=TRUE,
outline.col = "white",
ggtheme = ggplot2::theme_gray,
colors = c("#6D9EC1", "white", "#E46726")
)
Se puede apreciar cómo la variable atemp y temp estan muy correlacionadas, de ambas, seleccionaremos solamente la columna temp y eliminaremos la columna atemp para una mejora en la prediccion. Con la variable target count, estan directamente relacionadas las columnas temp, year, hour y windspeed y la columna humidity esta indirectamente realcionada.
Eliminamos la columna atemp del dataset.
p1_train <- p1_train %>% select(-atemp)
p1_test <- p1_test %>% select(-atemp)
predictors <- predictors %>% select(-atemp)
Calculamos la correlación no lineal MIC de las columnas con la variable target.
numpred <- predictors %>% select(where(is.numeric))
corrs <- data.frame(MIC = mine(numpred, y = p1_train$count, alpha = 0.7)) %>% select(MIC.Y)
corrs$Pearson <- cor(numpred, y = p1_train$count, method = 'pearson')
rownames(corrs) <- rownames(corrs$Pearson)
corrs
Podremos apreciar que las correlaciones no lineadas calculadas son (en proporcion) similares a las correlaciones de Pearson.
cbind(corrs$MIC.Y, corrs$MIC.Y ) %>% heatmaply_cor(xlab = "",
ylab = "Columnas", main = "Correlacion MIC de las columnas numericas con la variable target 'count'", cellnote = cbind(corrs$MIC.Y, corrs$MIC.Y ), k_col = 1, k_row = 1)
NA
El gràfico superior tiene 2 columnas pero se debe considerar una sola (estan repetidas), ha sido la unica manera que hemos encontrado de hacer el gráfico del MIC.
Realizamos un PCA sobre los datos para extraer los componentes principales y estudiarlos.
pca <- predictors %>% select(where(is.numeric)) %>% sample_n(200) %>%
prcomp(scale. = TRUE, center = TRUE)
summary(pca)
Importance of components:
PC1 PC2 PC3 PC4 PC5
Standard deviation 1.2906 1.0555 0.9785 0.8406 0.7458
Proportion of Variance 0.3331 0.2228 0.1915 0.1413 0.1112
Cumulative Proportion 0.3331 0.5560 0.7474 0.8888 1.0000
pca %>%
ggbiplot::ggbiplot(scale = 1)
Gràfico que representa la ponderación de las variables utilizadas para el PCA.
La variable target ‘count’ es una variable numerica que representa el numero de bicicletas alquiladas en las grandes ciudades. La variable tiene una media de 190 y un std de 182, un valor bastante elevado que indica mucha variación en la variable.
corrs$Pearson
[,1]
year 0.2623084
hour 0.4037404
temp 0.3971578
humidity -0.3245222
windspeed 0.1078587
Correlaciones de pearson de las columnas predictivas con la variable ‘count’.
mean(p1_train$count)
[1] 191.5154
sd(p1_train$count)
[1] 182.139
El grafico siguiente es un boxplot y un violin plot sobre la variable target ‘count’. Por el violin plot podemos ver que la mayoria de valores se encuentran cercanos al 0 pero por el boxplot vemos que los valores altos del dataset hacen que la media suba considerablemente. En el boxplot podemos apreciar también diversos outliers que sobresalen a partir del valor 550.
p1_train %>%
ggplot( aes(x=count, y=count)) +
geom_violin(width=1, color = "gray", alpha = 0.2, fill = 'green' ) +
geom_boxplot(width=0.3, color="black", fill='blue', alpha=0.3, outlier.colour="red", outlier.shape=8,
outlier.size=1, notch=TRUE) +
theme(
legend.position="none",
plot.title = element_text(size=11)
) +
ggtitle("Violin y box plot de los valores de la variable target ") +
xlab("") + theme_classic()
Partiendo de la base que la variable target es continua, podemos apreciar que la frequencia de la misma sigue una distribución parecida a la exponencial, ya que cómo más se aleja del número 0 menos frequencia encontramos. Habrá que tener en cuenta entonces, que el modelo va a tender a predecir valores bajos en la mayoria de los casos.
p1_train %>%
ggplot( aes(x=count)) +
geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth = 20)+
geom_density(alpha=.2, fill="#FF6666") + theme_classic()
Hacemos un escalado sobre los datos para normalizar.
p1_train <- p1_train %>% mutate(across(where(is.numeric), scale))
Entrenamos un modelo lineal
mod_lm <- lm(count~., p1_train)
Calculamos el mse sobre el train
mean(mod_lm$residuals^2) #mse
[1] 19913.19
Calculamos el msle sobre el train
mean(log(mod_lm$residuals^2)) #mse
[1] 8.389741
Se puede apreciar que el mse y msle es muy grande, el modelo no funciona bien.
plot(p1_train$count)
abline(mod_lm)
Se puede ver cómo un modelo lineal es demasiado sencillo para entender la complejidad de los datos y senzillamente predice la media del resultado.
Podremos coger el error 20000 (que representa el mse entre la media y cada valor) para poder decidir en los modelos posteriores si estan aprendiendo.
mod_glm <- glm(count~., p1_train, family = poisson())
mean(mod_glm$residuals^2)
[1] 0.7910352
Se ha intentado analizar la predicción pero no se encuentra sentido en la misma. El mse de 0.79 parece demasiado bajo para ser correcto y el predict sobre el train da valores count lejanos al ground truth.
set.seed(1234) #definimos una random seed para obtener los mismos resultados consistentemente
Definimos el modelo. Durante la realización de la práctica se ha hecho un grid search y se definen los parametros que devuelven mejor resultado.
rf <- randomForest(count~., data = p1_train, mtry = 10, importance = TRUE, ntree = 50 )
rf
Call:
randomForest(formula = count ~ ., data = p1_train, mtry = 10, importance = TRUE, ntree = 50)
Type of random forest: regression
Number of trees: 50
No. of variables tried at each split: 9
Mean of squared residuals: 2234.326
% Var explained: 93.26
Visualizamos que predictores son más importantes.
importance(rf)
%IncMSE IncNodePurity
year 111.315166 22945095.4
hour 186.937499 154124280.2
season 40.132400 9248360.1
holiday 7.192508 699533.4
workingday 103.924054 19701681.3
weather 19.992407 3985389.6
temp 43.378082 30699917.8
humidity 27.163578 8386444.8
windspeed 6.956011 3337028.5
varImpPlot(rf)
Podemos apreciar cómo la variable ‘hour’ es con diferencia la que mejor ayuda a predecir el numero de bicicletas alquiladas.
Analisis del model performance (mse) sobre el conjunto de train
mseDF <- data.frame(pred = rf$predicted, gt = rf$y)
mseDF
Podemos ver las diferencias entre el count predecido y el ground truth.
mean((mseDF$pred-mseDF$gt)^2) #mse
[1] 2234.326
Calculamos RMSLE (error cuadratico medio logaritmico)
rmsle(mseDF$gt,mseDF$pred)
[1] 0.3695941
Obtenemos un mse y rmsle coherentes, el modelo funciona.
Volvemos a hacer random forest con variables más correlacionadas con la columna a predecir.
p1_trainImp <- p1_train %>% select('hour','year','workingday','season','humidity', 'count')
Entrenamos random forest
rf2 <- randomForest(count~., data = p1_trainImp, mtry = 10, importance = TRUE, ntree = 50 )
rf2
Call:
randomForest(formula = count ~ ., data = p1_trainImp, mtry = 10, importance = TRUE, ntree = 50)
Type of random forest: regression
Number of trees: 50
No. of variables tried at each split: 5
Mean of squared residuals: 4122.14
% Var explained: 87.57
Comparación predicted vs ground truth.
mseDF2 <- data.frame(pred = rf2$predicted, gt = rf2$y)
mseDF2
Calculamos error
mean((mseDF2$pred-mseDF2$gt)^2) #mse
[1] 4122.14
Fucniona peor cogiendo solo las variables más importantes (mse 4122 > mse 2234). Nos quedaremos con la variable ’RF
Realizaremos una predicción sobre el dataset de test con el mejor modelo, en nuestro caso el RF.
p1_test <- p1_test %>% mutate(workingday = as.factor(workingday),
holiday = as.factor(holiday),
weather = as.factor(weather),
season = as.factor(season)) %>% mutate(across(where(is.numeric), scale))
Se realiza la prediccion sobre los datos de test con el modelo RandomForest.
predsFin <- predict(rf, p1_test)
data.frame(prediction_test = predsFin)
El algoritmo de regresion random forest es el que mejor predice el numero de bicicletas que se alquilaran un dia, con un mse de 2000 sobre el conjunto de train. Podemos ver analizando los resultados que el modelo no tiene mucho overfitting y generaliza bien, esto es asi en parte al haber creado un arbol no demasiado grande, evitando así la caída de valores en hojas.
Teniendo en cuenta que el mse entre la media y todos los valores del dataset era de 20.000 y con el mejor modelo hemos obtenido un mse de 2000, podemos concluir que nuestro modelo ha aprendido satisfactoriamente.
El rmsle definitivo es 0.3695941.
Cómo algoritmos de regresion se podrian intentar también el xgboost o el catboost pero perderiamos entendimiento de lo que esta haciendo el modelo, asi que se ha obviado para este trabajo.